6 Dielectric Properties

In this notebook, we will investigate the varying backscatter values associated with different land surfaces like water bodies, forests, grasslands and urban areas. We will use backscatter data from the Sentinel-1 satellite and we will utilize the CORINE Land Cover dataset to classify and extrapolate these surfaces, enabling us to analyze how different land cover types influence backscatter responses.

import json

import holoviews as hv
import intake
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from matplotlib.colors import BoundaryNorm, ListedColormap

from mrs.catalog import CorineColorCollection, get_intake_url
from mrs.plot import plot_corine_data, plot_variability_over_time

hv.extension("bokeh")

ⓘ

/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.11/abc.py:106: FutureWarning: xarray subclass DatasetHasLandCoverSig0 should explicitly define __slots__
  cls = super().__new__(mcls, name, bases, namespace, **kwargs)
/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.11/abc.py:106: FutureWarning: xarray subclass DatasetHasIntensity should explicitly define __slots__
  cls = super().__new__(mcls, name, bases, namespace, **kwargs)
/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.11/abc.py:106: FutureWarning: xarray subclass DatasetHasIntensityPhase should explicitly define __slots__
  cls = super().__new__(mcls, name, bases, namespace, **kwargs)
/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.11/abc.py:106: FutureWarning: xarray subclass DatasetHasBandData should explicitly define __slots__
  cls = super().__new__(mcls, name, bases, namespace, **kwargs)
/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.11/abc.py:106: FutureWarning: xarray subclass DatasetHasTopoPhase should explicitly define __slots__
  cls = super().__new__(mcls, name, bases, namespace, **kwargs)

6.1 Load Sentinel-1 Data

For our analysis we are using sigma naught backscatering data from Sentinel-1. The images we are analyzing cover the region south of Vienna and west of Lake Neusiedl. We load the data and and apply again a preprocessing function. Here we extract the scaling factor and the date the image was taken from the metadata. We will focus our attention to a smaller area containing a part of the Lake Neusiedl Lake and its surrounding land. The obtainedxarray dataset and is then converted to an array, because we only have one variable, the VV backscatter values.

url = get_intake_url()
cat = intake.open_catalog(url)
sig0_da = cat.neusiedler.read().sig0.compute()

https://git.geo.tuwien.ac.at/public_projects/microwave-remote-sensing/-/raw/main/microwave-remote-sensing.yml

Let’s have a look at the data by plotting the first timeslice.

sig0_da.isel(time=0).plot(robust=True, cmap="Greys_r").axes.set_aspect("equal")

6.2 Load CORINE Landcover Data

We will load the CORINE Land Cover, which is a pan-European land cover and land use inventory with 44 thematic classes. The resolution of this classification is 100 by 100m and the file was created in 2018 (CORINE Land Cover).

cor_da = cat.corine.read().land_cover.compute()

6.2.1 Colormapping and Encoding

For the different land cover types we use the official color encoding which can be found in CORINE Land Cover.

# Load encoding and validate with Pydantic
with cat.corine_cmap.read()[0] as f:
    corine_colors = CorineColorCollection.model_validate(json.load(f))

# Create cmap and norm for plotting
colors = [info["color"].as_hex() for info in corine_colors.items]
categories = [info["value"] for info in corine_colors.items]
cmap = ListedColormap(colors)
norm = BoundaryNorm([*categories, max(categories) + 1], len(categories))

Now we can plot the CORINE Land Cover dataset.

present_landcover_codes = np.unique(cor_da.values[~np.isnan(cor_da.values)].astype(int))

# Create the plot
plot_corine_data(cor_da, cmap, norm, corine_colors, present_landcover_codes)

Now we are ready to merge the backscatter data (sig0_da) with the land cover dataset (cor_da) to have one dataset combining all data.

var_ds = xr.merge([sig0_da, cor_da]).drop_vars("band")
var_ds

/tmp/ipykernel_2719/125755716.py:1: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  var_ds = xr.merge([sig0_da, cor_da]).drop_vars("band")

<xarray.Dataset> Size: 108MB
Dimensions:      (time: 8, x: 1230, y: 1221)
Coordinates:
  * time         (time) datetime64[ns] 64B 2023-08-05T16:51:22 ... 2023-10-28...
  * x            (x) float64 10kB 5.282e+06 5.282e+06 ... 5.294e+06 5.294e+06
  * y            (y) float64 10kB 1.571e+06 1.571e+06 ... 1.559e+06 1.559e+06
    spatial_ref  int64 8B 0
Data variables:
    sig0         (time, y, x) float64 96MB -6.99 -7.32 -8.78 ... -14.32 -14.22
    land_cover   (y, x) float64 12MB 12.0 12.0 12.0 12.0 ... 41.0 41.0 41.0 41.0

xarray.Dataset

• Dimensions:
  • time: 8
  • x: 1230
  • y: 1221
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2023-08-05T16:51:22 ... 2023-10-...

    array(['2023-08-05T16:51:22.000000000', '2023-08-17T16:51:22.000000000',
           '2023-08-29T16:51:23.000000000', '2023-09-10T16:51:24.000000000',
           '2023-09-22T16:51:24.000000000', '2023-10-04T16:51:24.000000000',
           '2023-10-16T16:51:24.000000000', '2023-10-28T16:51:24.000000000'],
          dtype='datetime64[ns]')

  • x
    (x)
    float64
    5.282e+06 5.282e+06 ... 5.294e+06

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([5281995., 5282005., 5282015., ..., 5294265., 5294275., 5294285.],
          shape=(1230,))

  • y
    (y)
    float64
    1.571e+06 1.571e+06 ... 1.559e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([1570955., 1570945., 1570935., ..., 1558775., 1558765., 1558755.],
          shape=(1221,))

  • spatial_ref
    ()
    int64
    0

    GeoTransform :

    5281990.0 10.0 0.0 1570960.0 0.0 -10.0

    crs_wkt :

    PROJCS["Azimuthal_Equidistant",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Azimuthal_Equidistant"],PARAMETER["latitude_of_center",53],PARAMETER["longitude_of_center",24],PARAMETER["false_easting",5837287.81977],PARAMETER["false_northing",2121415.69617],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    false_easting :

    5837287.81977

    false_northing :

    2121415.69617

    geographic_crs_name :

    WGS 84

    grid_mapping_name :

    azimuthal_equidistant

    horizontal_datum_name :

    World Geodetic System 1984

    inverse_flattening :

    298.257223563

    latitude_of_projection_origin :

    53.0

    longitude_of_prime_meridian :

    0.0

    longitude_of_projection_origin :

    24.0

    prime_meridian_name :

    Greenwich

    projected_crs_name :

    Azimuthal_Equidistant

    reference_ellipsoid_name :

    WGS 84

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    spatial_ref :

    PROJCS["Azimuthal_Equidistant",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Azimuthal_Equidistant"],PARAMETER["latitude_of_center",53],PARAMETER["longitude_of_center",24],PARAMETER["false_easting",5837287.81977],PARAMETER["false_northing",2121415.69617],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    array(0)

• Data variables: (2)
  • sig0
    (time, y, x)
    float64
    -6.99 -7.32 -8.78 ... -14.32 -14.22

    array([[[ -6.99,  -7.32,  -8.78, ..., -26.06, -21.29, -18.64],
            [ -6.37,  -8.79,  -8.78, ..., -26.2 , -21.8 , -21.18],
            [ -6.96,  -8.36,  -9.  , ..., -26.2 , -24.44, -24.31],
            ...,
            [ -7.82,  -4.94,  -4.8 , ...,  -8.78,  -8.02,  -7.3 ],
            [ -7.41,  -9.33,  -8.67, ..., -10.69,  -6.36,  -5.22],
            [-11.11, -11.8 , -11.55, ..., -14.89,  -6.95,  -3.89]],
    
           [[ -8.5 , -10.65, -12.07, ..., -22.27, -20.91, -21.92],
            [ -9.  , -10.03, -10.82, ..., -16.46, -16.9 , -21.84],
            [ -6.74,  -7.87,  -8.1 , ..., -17.41, -22.64, -23.01],
            ...,
            [ -8.74,  -7.44,  -6.69, ..., -10.85,  -8.39,  -7.08],
            [ -9.39,  -9.55,  -8.9 , ..., -13.  , -11.03,  -9.15],
            [-10.86, -10.84,  -9.7 , ..., -17.47, -14.58, -10.65]],
    
           [[-10.39, -11.92, -11.74, ..., -19.02, -20.28, -17.92],
            [ -8.65, -14.64, -12.34, ..., -21.77, -19.55, -16.67],
            [ -8.11, -11.29, -10.25, ..., -23.16, -20.23, -16.73],
            ...,
    ...
            [ -9.82,  -9.11,  -8.73, ...,  -9.79,  -7.33,  -6.92],
            [-12.23, -11.64,  -9.59, ..., -12.29,  -7.65,  -8.14],
            [-13.6 , -13.41, -11.73, ..., -19.82, -12.53, -11.83]],
    
           [[-12.49, -12.55, -12.89, ..., -21.35, -23.29, -23.87],
            [-12.15, -12.38, -14.03, ..., -18.07, -23.53, -23.23],
            [-12.41, -12.52, -14.08, ..., -18.7 , -21.47, -22.17],
            ...,
            [ -6.87,  -5.98,  -5.32, ..., -18.47, -12.26, -11.2 ],
            [-11.18,  -6.31,  -6.12, ..., -20.17, -12.56, -11.61],
            [-11.98,  -9.02,  -8.46, ..., -13.06, -10.4 ,  -9.  ]],
    
           [[ -6.41,  -6.59,  -7.16, ..., -22.25, -25.96, -25.45],
            [ -6.75,  -9.42,  -9.54, ..., -20.65, -26.34, -23.53],
            [ -8.24, -12.04, -10.46, ..., -20.65, -20.63, -21.05],
            ...,
            [ -6.64,  -5.06,  -4.89, ..., -13.67, -12.26,  -9.96],
            [ -9.62,  -8.06, -10.02, ..., -16.11, -14.74, -10.64],
            [-11.93, -11.41, -12.68, ..., -19.01, -14.32, -14.22]]],
          shape=(8, 1221, 1230))

  • land_cover
    (y, x)
    float64
    12.0 12.0 12.0 ... 41.0 41.0 41.0

    AREA_OR_POINT :

    Area

    DataType :

    Thematic

    RepresentationType :

    THEMATIC

    STATISTICS_COVARIANCES :

    136.429646247598

    STATISTICS_MAXIMUM :

    48

    STATISTICS_MEAN :

    25.753373398066

    STATISTICS_MINIMUM :

    1

    STATISTICS_SKIPFACTORX :

    1

    STATISTICS_SKIPFACTORY :

    1

    STATISTICS_STDDEV :

    11.680310194836

    array([[12., 12., 12., ..., 41., 41., 41.],
           [12., 12., 12., ..., 41., 41., 41.],
           [12., 12., 12., ..., 41., 41., 41.],
           ...,
           [ 2.,  2.,  2., ..., 41., 41., 41.],
           [ 2.,  2.,  2., ..., 41., 41., 41.],
           [ 2.,  2.,  2., ..., 41., 41., 41.]], shape=(1221, 1230))

6.3 Backscatter Variability

With this combined dataset we can study backscatter variability in relation to natural media. For example we can look at the backscatter variability for water by clipping the dataset to only contain the land cover class water, like so:

# 41 = encoded value for water bodies
encoded_value_for_waterbodies = 41
waterbodies_mask = var_ds.land_cover == encoded_value_for_waterbodies
waterbodies_mask.plot().axes.set_aspect("equal")

This gives use backscatter values over water only.

waterbodies_sig0 = var_ds.sig0.isel(time=0).where(waterbodies_mask)
waterbodies_sig0.plot(robust=True, cmap="Greys_r").axes.set_aspect("equal")

To get an idea of the variability we can create a histogram. Radar backscatter from water bodies fluctuates with surface roughness, which changes with wind conditions, creating spatial and temporal variations in signal intensity.

waterbodies_sig0.plot.hist(bins=50, edgecolor="black")

(array([5.0000e+00, 5.0000e+00, 8.0000e+00, 1.4000e+01, 3.5000e+01,
        7.8000e+01, 1.5600e+02, 3.0500e+02, 5.4000e+02, 1.0380e+03,
        1.9550e+03, 3.6850e+03, 6.9910e+03, 1.2269e+04, 1.9088e+04,
        2.6859e+04, 3.6217e+04, 4.5276e+04, 5.0939e+04, 4.9020e+04,
        3.9729e+04, 2.5400e+04, 1.2426e+04, 5.4390e+03, 2.6990e+03,
        2.0500e+03, 2.1330e+03, 2.6890e+03, 3.1490e+03, 3.9510e+03,
        4.2970e+03, 4.1110e+03, 3.6350e+03, 2.5900e+03, 1.6020e+03,
        7.5000e+02, 3.5200e+02, 1.2300e+02, 5.6000e+01, 3.7000e+01,
        2.8000e+01, 2.1000e+01, 1.3000e+01, 9.0000e+00, 1.0000e+01,
        8.0000e+00, 1.0000e+01, 5.0000e+00, 6.0000e+00, 4.0000e+00]),
 array([-37.91  , -36.9202, -35.9304, -34.9406, -33.9508, -32.961 ,
        -31.9712, -30.9814, -29.9916, -29.0018, -28.012 , -27.0222,
        -26.0324, -25.0426, -24.0528, -23.063 , -22.0732, -21.0834,
        -20.0936, -19.1038, -18.114 , -17.1242, -16.1344, -15.1446,
        -14.1548, -13.165 , -12.1752, -11.1854, -10.1956,  -9.2058,
         -8.216 ,  -7.2262,  -6.2364,  -5.2466,  -4.2568,  -3.267 ,
         -2.2772,  -1.2874,  -0.2976,   0.6922,   1.682 ,   2.6718,
          3.6616,   4.6514,   5.6412,   6.631 ,   7.6208,   8.6106,
          9.6004,  10.5902,  11.58  ]),
 <BarContainer object of 50 artists>)

6.4 Variability over Time

Next we will look at the changes in variability in backscatter values over time for each of the CORINE Land Cover types. We do this by creating the following interactive plot. We can spot that backscatter in agricultural fields varies due to seasonal cycles like planting, growing, and harvesting, each of which changes vegetation structure. Changes in backscatter are strongly related to soil moisture content from irrigation or rainfall. Ultimately, phenological stages of crops and canopy moisture dynamics can affect the backscatter signal.

plot_variability_over_time(corine_colors.to_dict(), var_ds, present_landcover_codes)